
# ==============================================================================
# DESCRIPTION ----
# ==============================================================================

# Conducts time-series cross sectional matching analysis of amakudari hires on 
# government loans received by private sector firms. 

# ==============================================================================
# LIBRARIES AND IMPORT ----
# ==============================================================================

# Options
options(scipen=999)

# Libraries
library(tidyverse)
library(DeclareDesign)
library(PanelMatch)

# Functions
source("code/0.2_functions.R")

# Data import
load("data/loans_amakudari.Rdata")

# ==============================================================================
# FINAL DATA PREP ----
# ==============================================================================
# Remove duplicates in data ----------------------------------------------------
financials <- loans_covs %>% 
  distinct(firm_dest_en, year, .keep_all = TRUE) %>%
  filter(year != 2020)

# Create alternative amakudari coding where 1 for all periods after hire -------
fin <- financials %>% 
  mutate(amakudari_fill = na_if(as.numeric(amakudari_binary), 0),
         amakudari_match = na_if(as.numeric(amakudari_binary), 0)) %>%
  group_by(nikkei_code) %>%
  fill(amakudari_fill, .direction = "down") %>%
  fill(amakudari_match, .direction = "downup") %>%
  ungroup() %>%
  mutate(amakudari_fill = as.numeric(replace_na(amakudari_fill, 0)),
         amakudari_match = as.numeric(replace_na(amakudari_match, 0)),
         year = as.integer(as.numeric((as.factor(year)))),
         nikkei_code = as.numeric(factor(nikkei_code)))
fin <- as.data.frame(fin)

# Create amakudari codings for each ministry -----------------------------------
ministry <- financials %>% 
  mutate(
    amakudari_binary_METI_MOF = amakudari_binary_METI + amakudari_binary_MOF,
    amakudari_binary_METI_MOF = ifelse(amakudari_binary_METI_MOF >= 1, 1, amakudari_binary_METI_MOF),
    amakudari_binary_other = amakudari_binary_MIAC + amakudari_binary_MLIT + 
      amakudari_binary_MHLW + amakudari_binary_MAFF + amakudari_binary_MEXT +
      amakudari_binary_MOD + amakudari_binary_MOE + amakudari_binary_MOFA +
      amakudari_binary_MOJ,
    amakudari_binary_other = ifelse(amakudari_binary_other >= 1, 1, amakudari_binary_other),
    across(starts_with("amakudari"), ~na_if(as.numeric(.), 0))
  ) %>%
  group_by(nikkei_code) %>%
  fill(c(amakudari_n_Other:amakudari_binary_other), .direction = "down") %>%
  ungroup() %>%
  mutate(
    across(starts_with("amakudari"), ~as.numeric(replace_na(., 0))),
    year = as.integer(as.numeric((as.factor(year)))),
    nikkei_code = as.numeric(factor(nikkei_code))
  )
ministry <- as.data.frame(ministry)

# ==============================================================================
# VISUALIZE DATA ----
# ==============================================================================
# Define panelmatch data
loan.panel <- PanelData(panel.data = fin,
                        unit.id = "nikkei_code",
                        time.id = "year",
                        treatment = "amakudari_fill",
                        outcome = "debt_total_public")

# FIGURE A9: Visualize treatment
cont_treat <- DisplayTreatment(panel.data = loan.panel,
  legend.position = "none", xlab = "Year", ylab = "Firm", y.size = NULL,
  hide.y.tick.label = TRUE, dense.plot = TRUE,
  color.of.treated = "seagreen3", color.of.untreated = "steelblue2")

ggsave("figures/a9_loan_treat_control.pdf", cont_treat, height = 5, width = 6)

# ==============================================================================
# PRIMARY ESTIMATES: TIME SERIES CROSS SECTIONAL MATCHING ----
# ==============================================================================
# All ministries ---------------------------------------------------------------
loan_result <- run_panelmatch(
  data = fin, treatment = "amakudari_fill", outcome = "debt_total_public")

# METI -------------------------------------------------------------------------
loan_result_meti <- run_panelmatch(
  data = ministry, treatment = "amakudari_binary_METI", outcome = "debt_total_public")

# MOF --------------------------------------------------------------------------
loan_result_mof <- run_panelmatch(
  data = ministry, treatment = "amakudari_binary_MOF", outcome = "debt_total_public")

# METI and MOF -----------------------------------------------------------------
loan_result_meti_mof <- run_panelmatch(
  data = ministry, treatment = "amakudari_binary_METI_MOF", outcome = "debt_total_public")

# ALL OTHER MINISTRIES ---------------------------------------------------------
loan_result_other <- run_panelmatch(
  data = ministry, treatment = "amakudari_binary_other", outcome = "debt_total_public")

# ==============================================================================
# ROBUSTNESS CHECKS: TIME SERIES CROSS SECTIONAL MATCHING ----
# ==============================================================================
# Expand to all covariates that are plausibly pre-treatment --------------------
loan_result_all <- run_panelmatch(
  data = fin, treatment = "amakudari_fill", outcome = "debt_total_public",
  covariates = ~ total_assets + total_liabilities + employees + ebitda + 
    gross_profit + operating_revenue + leverage + reserve_ratio + roi + roe)

# Examine effects on private loans --------------------------------------------
loan_result_private <- run_panelmatch(
  data = fin, treatment = "amakudari_fill", outcome = "debt_total_private")

# 2 lag periods (all ministries) -----------------------------------------------
loan_result_lag <- run_panelmatch(
  data = fin, treatment = "amakudari_fill", outcome = "debt_total_public", 
  lag = 2)

# 2 lag periods (METI AND MOF) -------------------------------------------------
loan_result_lag_meti_mof <- run_panelmatch(
  data = ministry, treatment = "amakudari_binary_METI_MOF", 
  outcome = "debt_total_public", lag = 2)

# Covariate Balancing Propensity Score matching: all ministries ----------------
loan_result_cbps <- run_panelmatch(
  data = fin, treatment = "amakudari_fill", outcome = "debt_total_public",
  refinement.method = "CBPS.match")

# Covariate Balancing Propensity Score matching: METI and MOF only -------------
loan_result_cbps_meti_mof <- run_panelmatch(
  data = ministry, treatment = "amakudari_binary_METI_MOF", 
  outcome = "debt_total_public", refinement.method = "CBPS.match")

# Covariate Balancing Propensity Score matching (MSM): all ministries ----------
loan_result_cbps_msm <- run_panelmatch(
  data = fin, treatment = "amakudari_fill", outcome = "debt_total_public",
  refinement.method = "CBPS.msm.weight")

# Covariate Balancing Propensity Score matching (MSM): METI and MOF only -------
loan_result_cbps_msm_meti_mof <- run_panelmatch(
  data = ministry, treatment = "amakudari_binary_METI_MOF", 
  outcome = "debt_total_public", refinement.method = "CBPS.msm.weight")

# Traditional propensity score matching: all ministries ------------------------
loan_result_ps <- run_panelmatch(
  data = fin, treatment = "amakudari_fill", outcome = "debt_total_public",
  refinement.method = "ps.match")

# Traditional propensity score matching: METI and MOF only ---------------------
loan_result_ps_meti_mof <- run_panelmatch(
  data = ministry, treatment = "amakudari_binary_METI_MOF", 
  outcome = "debt_total_public", refinement.method = "ps.match")

# Different numbers of lead periods --------------------------------------------
# 3 leads
loan_result_3lead <- run_panelmatch(data = fin, 
  treatment = "amakudari_fill", outcome = "debt_total_public", lead = 0:3)

# 4 leads
loan_result_4lead <- run_panelmatch(data = fin, 
  treatment = "amakudari_fill", outcome = "debt_total_public", lead = 0:4)

# 6 leads
loan_result_6lead <- run_panelmatch(data = fin, 
  treatment = "amakudari_fill", outcome = "debt_total_public", lead = 0:6)

# 7 leads
loan_result_7lead <- run_panelmatch(data = fin, 
  treatment = "amakudari_fill", outcome = "debt_total_public", lead = 0:7)

# ==============================================================================
# OUTPUT FIGURES AND TABLES: PRIMARY ESTIMATES ----
# ==============================================================================
# All ministries ---------------------------------------------------------------
# FIGURE 5A: Plot results and save
ggplot(loan_result$cleaned_result_table) + gglayers_update +
  scale_y_continuous(limits = c(-5000, 20000), breaks= seq(-5000, 20000, 5000))
ggsave("figures/5a_tscs_loan.pdf", width = 10, height = 5)

# TABLE A6: Output table
stargazer::stargazer(
  loan_result$cleaned_result_table, 
  summary = FALSE, digits = 2, rownames = FALSE,
  title = "Estimated effect of bureaucratic hires on size of government loans received, by year afer hire",
  label = "tab: tscs_loan",
  out = "tables/a6_tscs_loan.tex",
  notes = paste("Note: Matched sets =", length(loan_result$estimate_object$matched.sets))
)

# METI ---------------------------------------------------------------
# FIGURE A10: Plot results and save
ggplot(loan_result_meti$cleaned_result_table) + gglayers_update +
  scale_y_continuous(limits = c(-10000, 25000), breaks= seq(-10000, 25000, 5000))
ggsave("figures/a10_tscs_loan_meti.pdf", width = 10, height = 5)

# TABLE A7: Output to LaTeX using stargazer
stargazer::stargazer(
  loan_result_meti$cleaned_result_table, 
  summary = FALSE, digits = 2, rownames = FALSE,
  title = "Estimated effect of bureaucratic hires on size of government loans received, METI re-hires only",
  label = "tab: tscs_loan_meti",
  out = "tables/a7_tscs_loan_meti.tex",
  notes = paste("Note: Matched sets =", 
                length(loan_result_meti$estimate_object$matched.sets))
)

# MOF --------------------------------------------------------------------------
# FIGURE A11: Plot results and save
ggplot(loan_result_mof$cleaned_result_table) + gglayers_update +
  scale_y_continuous(limits = c(-10000, 25000), breaks= seq(-10000, 25000, 5000))
ggsave("figures/a11_tscs_loan_mof.pdf", width = 10, height = 5)

# TABLE A8: Output to LaTeX using stargazer
stargazer::stargazer(
  loan_result_mof$cleaned_result_table, summary = FALSE, digits = 2, rownames = FALSE,
  title = "Estimated effect of bureaucratic hires on size of government loans received, MOF re-hires only",
  label = "tab: tscs_loan_mof",
  out = "tables/a8_tscs_loan_mof.tex",
  notes = paste("Note: Matched sets =", length(loan_result_mof$estimate_object$matched.sets))
)

# METI AND MOF -----------------------------------------------------------------
# FIGURE 5B: Plot results and save
ggplot(loan_result_meti_mof$cleaned_result_table) + gglayers_update +
  scale_y_continuous(limits = c(-5000, 20000), breaks= seq(-5000, 20000, 5000))
ggsave("figures/5b_tscs_loan_meti_mof.pdf", width = 10, height = 5)

# TABLE A9: Output to LaTeX using stargazer
stargazer::stargazer(
  loan_result_meti_mof$cleaned_result_table, 
  summary = FALSE, digits = 2, rownames = FALSE,
  title = "Estimated effect of bureaucratic hires on size of government loans received, METI and MOF re-hires only",
  label = "tab: tscs_loan_meti_mof",
  out = "tables/a9_tscs_loan_meti_mof.tex",
  notes = paste("Note: Matched sets =", length(loan_result_meti_mof$estimate_object$matched.sets))
)

# All ministries other than METI and MOF ---------------------------------------
# FIGURE A12: Plot results and save
ggplot(loan_result_other$cleaned_result_table) + gglayers_update +
  scale_y_continuous(limits = c(-15000, 10000), breaks= seq(-15000, 10000, 5000))
ggsave("figures/a12_tscs_loan_other.pdf", width = 10, height = 5)

# TABLE A10: Output to LaTeX using stargazer
stargazer::stargazer(
  loan_result_other$cleaned_result_table, 
  summary = FALSE, digits = 2, rownames = FALSE,
  title = "Estimated effect of bureaucratic hires on size of government loans received, all ministries other than METI and MOF",
  label = "tab: tscs_loan_other",
  out = "tables/a10_tscs_loan_other.tex",
  notes = paste("Note: Matched sets =", length(loan_result_other$estimate_object$matched.sets))
)

# ==============================================================================
# OUTPUT FIGURES AND TABLES: ROBUSTNESS CHECKS ----
# ==============================================================================
# FIGURE A18: All covariates (including plausibly post-treatment) --------------
ggplot(loan_result_all$cleaned_result_table) + gglayers_update
ggsave("figures/a18_tscs_loan_posttreatment.pdf", width = 10, height = 5)

# Private loans ----------------------------------------------------------------
# FIGURE A19: Plot results and save
ggplot(loan_result_private$cleaned_result_table) + gglayers_update + 
  scale_y_continuous(limits = c(-35000, 20000), 
                     breaks= seq(-35000, 20000, 5000))
ggsave("figures/a19_tscs_loan_private.pdf", width = 10, height = 5)

# TABLE A13: Output to LaTeX using stargazer
stargazer::stargazer(
  loan_result_private$cleaned_result_table, summary = FALSE, digits = 2, rownames = FALSE,
  title = "Estimated effect of bureaucratic hires on size of private loans received",
  label = "tab: tscs_loan_private",
  out = "tables/a13_tscs_loan_private.tex",
  notes = paste("Note: Matched sets =", length(loan_result_private$estimate_object$matched.sets))
)

# 2 lag periods (all ministries) -----------------------------------------------
# FIGURE 14: Plot results and save
ggplot(loan_result_lag$cleaned_result_table) + gglayers_update +
  scale_y_continuous(limits = c(-5000, 15000), breaks= seq(-5000, 15000, 5000))
ggsave("figures/a14_tscs_loan_2lags.pdf", width = 10, height = 5)

# TABLE 11: Output to LaTeX using stargazer
stargazer::stargazer(
  loan_result_lag$cleaned_result_table, 
  summary = FALSE, digits = 2, rownames = FALSE,
  title = "Estimated effect of bureaucratic hires on size of government loans received, requiring 2 lag periods",
  label = "tab: tscs_loan_lag",
  out = "tables/a11_tscs_loan_lag.tex",
  notes = paste("Note: Matched sets =", length(loan_result_lag$estimate_object$matched.sets))
)

# 2 lag periods (METI AND MOF) -------------------------------------------------
# FIGURE A16B: Plot results and save
ggplot(loan_result_lag_meti_mof$cleaned_result_table) + gglayers_update +
  scale_y_continuous(limits = c(-5000, 15000), breaks= seq(-5000, 15000, 5000))
ggsave("figures/a16b_tscs_loan_2lags_meti_mof.pdf", width = 10, height = 5)

# TABLE A12: Output to LaTeX using stargazer
stargazer::stargazer(
  loan_result_lag_meti_mof$cleaned_result_table, summary = FALSE, digits = 2, rownames = FALSE,
  title = "Estimated effect of bureaucratic hires on size of government loans received, METI and MOF re-hires only, requiring 2 lag periods",
  label = "tab: tscs_loan_lag_meti_mof",
  out = "tables/a12_tscs_loan_lag_meti_mof.tex",
  notes = paste("Note: Matched sets =", length(loan_result_lag_meti_mof$estimate_object$matched.sets))
)

# Covariate Balancing Propensity Score matching: all ministries ----------------
# FIGURE A20C: Plot results and save
ggplot(loan_result_cbps$cleaned_result_table) + gglayers_update +
  scale_y_continuous(limits = c(-5000, 15000), breaks= seq(-5000, 15000, 5000))
ggsave("figures/a20c_tscs_loan_cbps.pdf", width = 10, height = 5)

# Covariate Balancing Propensity Score matching: METI and MOF only -------------
# FIGURE A21C: Plot results and save
ggplot(loan_result_cbps_meti_mof$cleaned_result_table) + gglayers_update +
  scale_y_continuous(limits = c(-10000, 20000), breaks= seq(-10000, 20000, 5000))

ggsave("figures/a21c_tscs_loan_cbps_meti_mof.pdf", width = 10, height = 5)

# Covariate Balancing Propensity Score matching (MSM): all ministries ----------
# FIGURE A20d
ggplot(loan_result_cbps_msm$cleaned_result_table) + gglayers_update +
  scale_y_continuous(limits = c(-5000, 15000), breaks= seq(-5000, 15000, 5000))
ggsave("figures/a20d_tscs_loan_cbps_msm.pdf", width = 10, height = 5)

# Covariate Balancing Propensity Score matching (MSM): METI and MOF only -------
# FIGURE A21d
ggplot(loan_result_cbps_msm_meti_mof$cleaned_result_table) + gglayers_update +
  scale_y_continuous(limits = c(-10000, 15000), breaks= seq(-10000, 15000, 5000))
ggsave("figures/a21d_tscs_loan_cbps_msm_meti_mof.pdf", width = 10, height = 5)

# Traditional propensity score matching: all ministries ------------------------
# FIGURE A20b
ggplot(loan_result_ps$cleaned_result_table) + gglayers_update +
  scale_y_continuous(limits = c(-5000, 15000), breaks= seq(-5000, 15000, 5000))
ggsave("figures/a20b_tscs_loan_ps.pdf", width = 10, height = 5)

# Traditional propensity score matching: METI and MOF only ---------------------
# FIGURE A21b
ggplot(loan_result_ps_meti_mof$cleaned_result_table) + gglayers_update +
  scale_y_continuous(limits = c(-5000, 15000), breaks= seq(-5000, 15000, 5000))
ggsave("figures/a21b_tscs_loan_ps_meti_mof.pdf", width = 10, height = 5)

# Different numbers of lead periods --------------------------------------------
# Plot each lead plot
lead3 <- ggplot(loan_result_3lead$cleaned_result_table) + gglayers_leads
lead4 <- ggplot(loan_result_4lead$cleaned_result_table) + gglayers_leads
lead6 <- ggplot(loan_result_6lead$cleaned_result_table) + gglayers_leads
lead7 <- ggplot(loan_result_7lead$cleaned_result_table) + gglayers_leads

# FIGURE A13: Group plots together and export
ggpubr::ggarrange(lead3, lead4, lead6, lead7, ncol = 2, nrow = 2)
ggsave("figures/a13_tscs_loan_leads.pdf", width = 10, height = 8)

# ==============================================================================
# REFINEMENT PLOTS ----
# ==============================================================================
# Calculate covariate balance by refinement method -----------------------------
# All ministries
covbal <- get_covariate_balance(
  loan_result$match_object, loan_result_cbps$match_object,
  loan_result_cbps_msm$match_object, loan_result_ps$match_object,
  panel.data = loan_result$panel_object,
  covariates = c("total_assets", "total_liabilities", "employees", "ebitda", 
                 "operating_revenue", "gross_profit"),
  include.unrefined = TRUE)

# METI and MOF only
covbal_meti_mof <- get_covariate_balance(
  loan_result_meti_mof$match_object, loan_result_cbps_meti_mof$match_object,
  loan_result_cbps_msm_meti_mof$match_object, loan_result_ps_meti_mof$match_object,
  panel.data = loan_result_meti_mof$panel_object,
  covariates = c("total_assets", "total_liabilities", "employees", "ebitda", 
                 "operating_revenue", "gross_profit"),
  include.unrefined = TRUE)

# Create refinement plots ------------------------------------------------------
# FIGURE A23: All ministries
pdf("figures/a23_tscs_balance.pdf", width = 8, height = 8)
par(mfrow = c(2, 2))
plot(covbal[1], type = "scatter", ylab = "After refinement",
     main = "Mahalanobis")
plot(covbal[2], type = "scatter", ylab = "After refinement",
     main = "Covariate balanced propensity score")
plot(covbal[3], type = "scatter", ylab = "After refinement",
     main = "Covariate balanced proensity score (MSM)")
plot(covbal[4], type = "scatter", ylab = "After refinement",
     main = "Propensity score")
dev.off()

# FIGURE A24: METI and MOF only
pdf("figures/a24_tscs_balance_meti_mof.pdf", width = 8, height = 8)
par(mfrow = c(2, 2))
plot(covbal_meti_mof[1], type = "scatter", ylab = "After refinement",
     main = "Mahalanobis")
plot(covbal_meti_mof[2], type = "scatter", ylab = "After refinement",
     main = "Covariate balanced propensity score", )
plot(covbal_meti_mof[3], type = "scatter", ylab = "After refinement",
     main = "Covariate balanced proensity score (MSM)")
plot(covbal_meti_mof[4], type = "scatter", ylab = "After refinement",
     main = "Propensity score")
dev.off()
